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ABSTRACT 


This dissertation develops a new general method of solving Prony’s problem. Two 
special cases of this new method have been developed previously. They are the 
Matrix Pencil and the Osculatory Interpolation. The dissertation shows that they 
are instances of a more general solution type which allows a wide ranging class of linear 
functional to be used in the solution of the problem. This class provides a continuum 
of functionals which provide new methods that can be used to solve Prony’s problem. 
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CHAPTER 1 
INTRODUCTION 


Many real-world processes obey constant coefficient differential equations. The 
measurable quantities from these processes are described by exponential functions. 
Certain of these processes, such as radioactive decay, demand to be modeled using 
real decay factors. 

The scientist is often presented with the task of identifying the material in question 
using whatever measurement capabilities are at his disposal. The rate of radiation 
from a mass is one of the measureable quantities available. This leads to the task 
of determining the real decay factors from measured data. The problem then, is to 
identify the real coefficients bk in a real-valued function of the form 

= (i-i) 

k= 1 

from data. In order to make the problem more tractable, the data are acquired at 
uniform time intervals. This problem was first considered and a solution presented 
by Baron R. deProny, a Parisian engineer, in 1795 [85]. 

Prony had limited success in solving the problem. This is because the problem is 
highly nonlinear. In addition, unlike the problem of identifying a function with sums 
of sines and cosines, the basis functions of this problem, the real exponentials, are 
not orthogonal [56]. This has resulted in repeated attempts at solving the problem 
since Prony [89, 56, 29, 38, 73, 2, 39, 9, 83], 

In recent years, this problem has generated considerable interest. Until recently, 
the solutions have fallen into two categories: nonlinear gradient descent methods, 
and autoregressive methods. The gradient descent methods are very costly and don’t 
guarantee global convergence. Thus, most interest has focussed on the autoregressive 
or AR model. The use of the AR model stems from the observation that functions of 
the form (1.1) are, in the case of integer time intervals, solutions to some difference 
equation with constant coefficients. Only recently has this type of equation been 
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called autoregressive. However, a full theory has grown up around the AR problem. 
Fundamental to the theory is the assumption that the process is stationary. Distur- 
bances, known and unknown, are allowed in the process. Various levels of success 
have been reported for the AR approach to the problem. 

Recently, Ammar et al. [2] and Hua and Sarkar [39] developed a matrix pencil 
method for determining the coefficients b k . The matrix pencil method performs well 
at findin g the coefficients when noise levels are low, but is not as tolerant of noise as 
many of the AR methods. It is, however, directly related to the AR method proposed 

by Prony [79]. 

The pencil method and the oscillatory method developed by Martin et al. [65] 
are both specialized examples of a new class of methods for solving this problem. 
This new class of solutions is termed the linear functional method. By applying a 
linear functional to the data, the data are modified in such a way that the inverse of 
the functional applied to the eigenvalues of a certain generalized eigenvalue problem 
yield the desired coefficients. For the functionals tried so far, the functional method 
doesn’t perform well with noisy data. 

1.1 Background 

A time varying quantity is measured in many physical problems, for instance the 
temperature of a heated object at room temperature, or the rate of radiation from a 
radioactive object. 

In the case of a radioactive object, the rate of radiation is determined by the 
mixture of materials. 

Assumption 1. The mixtures of materials are such that the radiation rates remain in 
the linear region. 

Each type of radioactive material loses mass at a rate which is different than that of 
others. One task of interest then, is the determination of the different rates of decay 
of matter in the object. 
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Assume that the rate of decay of a homogenous mass is related to the amount of 
radioactive material in the object. Designate the rate of decay at a particular time, t, 
by p(t), the mass of remaining radioactive material by M(t) and let 7 be a constant, 

then 

p{t) = 7 M{t), 

and 

p(t) = M(i), 


so 

M(t) = 7 M{t). 


Thus 


M(t) = M 0 e 7t , 


(1.2) 


where A / 0 is the mass of the radioactive material at time t = 0 . 

Now, if the radioactive object is composed of more than one type of radioactive 
material, equation (1.2) is no longer theoretically sound. Assuming that the radiation 
from one material does not affect the rate of radiation from another material, the rates 
will be additive and the mass of the object at time t becomes 

M{t) = ^2a k e bkt , (!- 3 ) 

where the afc’s represent the initial mass of each of the component materials, and 
the b k ’s represent the decay rates for each of the materials in the mix. Of course the 
assumption that the mixture of materials does not affect the decay rates of the matter 
in the mix is false. However, in the mix, new decay rates will exist for each type of 
matter and if the constants are adjusted, the formula (1.3) is still representative of 
the behaviour of the object. By differentiation, 

p(t) = Y a* k e bkt , 

k= 1 

so that the problem is the same whether the mass or radiation is being measured. In 
both cases, the desired quantities are the b k ’s. 
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1.2 Roadmap 

The rest of the dissertation is organized as follows. The second chapter presents 
some fundamental results. These results establish the basic theory and notation which 
is used throughout the dissertation. Chapter 3 presents the linear functional method. 
This method is the primary result of the dissertation. Once it is in place, Chapter 
4 gives details on some possibilities for using the functional method. Chapter 5 
provides some numerical comparisons between some methods based on the functional 
technique and some standard methods which appear elsewhere. 
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CHAPTER 2 
FOUNDATION 


This chapter presents some basic results and prepares the groundwork for the 
rest of the dissertation. Some notational shortcuts are presented to make the later 
developments shorter. 


2.1 The Vandermonde 


Lemma 1. If Z\, z 2 , z$, . . . ,z p are distinct complex numbers, the Vandermonde matrix 


V = 


z 1 


Z2 


1 


2 3 

,2 

~3 


1 


■vP- 1 „P-1 „P- 1 


rP- 1 


is nonsingular. 


Proof. This standard result can be found in Atkinson [4], 

Lemma 2. If z\, z 2 , z$, . . . , z v are distinct and n > p, the Vandermonde matrix 

1 1 1 ... 1 


V = 


Z\ z 2 Z 3 


? 9 ? 

Zf 2 2 - Z 3 - 


z l z 2 z 3 * * ■ z p 


has full column rank. 

Proof. This follows from the previous lemma since V has rank p. ■ 

This dissertation focuses on numbers that can be represented by exponentials. 
Since any complex number except zero can be represented as an exponential, the zf s 
in the previous two lemmas are restricted to be nonzero complex numbers. Under 
this restriction, the lemmas still hold. 
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2.2 The Difference Equation 


Consider the function 


f(t) = a k e bkt . 
Jfc=l 


( 2 . 1 ) 


Assumption 2. The b k ’$ are distinct. 

Lemma 3. The functions f(t) of the form (2.1) satisfy a difference equation of the 
form 


f(t + n) + ctp-\f(t + n — 1) 4- • • • + cxof (t + n — p) — 0. (2-2) 

Furthermore, the roots, \ k) of the polynomial 

X p + oip-i X p 1 + ■ • • + q;o = g{ A) (2-3) 

are the exponential factors in (2.1). In other words, 

At = e bk . 


Proof. At time t + m, 

v 

f(t + m) = a k e bkt e bkTn (2.4) 

Jk=l 

and (2.1) becomes 





1 

P 

t— 

1 

/(() = [! 1 ••• 1 ] 

e bit 




e bpt 


i 

■£ 

i 


For f(t),f(t + l),...,f(t+p-l), 


m 


1 1 

1 


~e blt 



a i 

f{t + 1) 

= 

e bl e b - 

e bp 


e M 



a 2 

i 

+ 

1 ' 


e *i(?-l) e h(p-l) 

e b P (p- 1) 



e bpt 


Clp 


f{t ) =VT (t)a. (2.6) 
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Note that V is a Vandermonde matrix, and, as a direct result of the assumption, 
nonsingular. Thus, the vector a can be expressed as 

Then, according to (2.4) with the notation introduced in (2.6), 

f{t + 1) = VT(t + l)o = VT(£ + 1)T {ty'V-'fit). 

But using the definition of T(t) as a diagonal matrix, 

T{t + l)r(i )- 1 = T(l) = 

SO 

f{t + l) = VT(l)V- 1 f{t). (2.7) 

Indeed, the characteristic equation of the matrix T(l) is given by 

s (A) = ] ? [(A-e 6 ‘) = A'’ + ^a i A*. (2.8) 

fc=l k=0 

Let 

A = VT(1)V -1 . 

Then 
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and has the same characteristic equation. Writing out (2.7): 


7(*+i) 


m 

+ 2) 

=A 


f(t + P)_ 


j(t+p- 1)_ 


Observe that 


T 

O 

1 0 .. 

1 

O 


1 

io 

i 

0 

0 1 .. 

0 


f(t + l) 

_-aS -a; ... 



f{t + P~ 1)_ 


the matrix A is in companion form; thus the factors, 


matrix .4 are equal to the coefficients, a k , of its characteristic polynomial. 


(2.9) 


a * k , in the 
Designate 


f{n) = y n , then 


Un + Oip-\y n -\ H t" OtoVn-p — 0, 


( 2 . 10 ) 


which is a difference equation. Thus, the functions (2.1) satisfy the difference equa- 
tions (2.10). In addition, from (2.9), the coefficients b k of the function (2.1) are related 
to the coefficients, a k , of the difference equation by the relation (2.8). That is, if A, 
is a root of the polynomial g( A), then 

Ai = e bk (2-11) 


for one of the b k s. 


2.3 The Prony Step 

The preceding proof has a few details worth mentioning. First, the step of the 
derivation in which (2.6) was solved for a and the result put into (2.7) is very common 
in the derivations of this thesis. It is called the “Prony Step” because it effectively 
transforms a nonlinear problem into a linear one. The second aspect of the proof worth 
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noticing is that a diagonalization of the companion matrix is obtained provided the 
eigenvalues are distinct. If they are not distinct, then the Vandermonde is singular 
and this diagonalization does not hold. 

Lemma 4. Let 

0 1 0 ... 0 

0 0 1 ... 0 

G = 

— Oio —Oil ■ ■ ■ —Olp - 1 

If the eigenvalues, Ai, A 2 , . . . , A p , of G are distinct , then G = VW~ l , where 



and the columns ofV are eigenvectors of G. 

Proof. Consider GV. Without loss of generality consider the first column of the 
result. 
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But, since A! is a root of the characteristic equation of G, 


— do — oqAi ct 2 ^i — ‘ ‘ ' &p— 1^1 — Aj. 

Thus 




1 


T— < 

1 


1 

0 1 0 ... 0 

0 0 1 ... o 


Ai 

A? 

= 

A? 

A? 

= Ai 

Ai 

A? 

-do -d! . . . -ot v - 1 


1 

1 

i 


1 

^ .. 

i 


1 

1 


Then the first column of GV is the first column of V'A. The proof for the other 
columns is similar. Thus 


GV = V'A. 

■ 

2.4 Notation 

With these results in place, the standard notation to be used throughout this 
dissertation is presented. 

2.4.1 Model 

The model for the analysis is 

/(t)=y>e‘“. (2.12) 

k~ 1 

The t are taken to occur at the integers 0, 1,2,3, — Denote the function values at 
the times t = n by 

Vn = f{n). 
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2.4.2 Functionals 

In the following, examples will be given in which different functionals are used. 
The functionals will typically be given by a rule. In order to shorten the presentation, 
whenever there is no ambiguity, the values which result from applying the functionals 
to the data will be represented by y. For example, if the functional T is applied to 
the data, 

Vn ~ 3~(yn)- 


2.4.3 Hankel Matrices 

The following is the template for solution of the inverse problem of determining 
the coefficients a* of the difference equation 


y n + OCp-iVn-i + • • • + OiQy n -p = 0. 


Let 



y o 

y\ ■ • 

• 2/p- 1 

y i 

1/2 • • 

2/p 

y P - 1 

y P •• 

• l/2p-2 


be the square Hankel matrix formed from the data. Then, according to the difference 
equation, if 

— QO 

-Cq 

a = 

&p — 1 


and 


y P = 


Up 

y P + 1 


\V2p-l 


li 



then 


H Q a = y v . 

Hq is always the designation for the Hankel matrix of the data. It will sometimes be 
convenient to use Ho for the corresponding nonsquare Hankel matrix when there are 
more than 2 p data values. When needed, Hi represents the Hankel matrix made of 
values other than the data. 

For example, the Hankel matrix of data values starting at the second value, yi, 
instead of yo is 

yi y 2 ■■■ y P 

„ 2/2 2/3 ••• 2/p+i 

= • 

Up Vp + 1 • * • 3/2p — l 

Now, consider the solution problem which appears repeatedly in this dissertation. 
For some matrix Hi, solve the equation 

H 0 G = Hi, (2.13) 

and determine the eigenvalues of the matrix G. To make this process a little clearer, 
consider the result of multiplying the equation on the right by an eigenvector of G. 

HqGx = HiX 
H 0 \x = HiX 

{H 0 X - Hi)x = 0, (2.14) 

which implies that x and A are eigenvector and eigenvalue of the generalized eigenvalue 
equation (2.14). 

Now suppose that equation (2.14) has been solved for the eigenvalue and eigen- 
vector. Then 

HqXx = HiX. 
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If Ho is nonsingular, then this equation has the solution 


Ax = H 0 l H\X. 


But, in this case, (2.13) gives 


G = H^H,. 


Thus, the lemma: 

Lemma 5. If the matrix H 0 is nonsingular, then (A, x) are an eigenpair of the gener 
alized eigenvalue equation 

{Ho A - Hi)x = 0 x ^ 0 
if and only if (A, x) are an eigenpair of the matrix 

G = Hq 1 H]_. 

With this representation, it is obvious that it will not be necessary to determine 
the matrix G in order to find its eigenvalues. Similar results can be proven for 
the overdetermined case [106], In that case the appropriate generalized eigenvalue 

equation is 

(HjH 0 \ - HjH^x = 0 x^O. 
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CHAPTER 3 

THE LINEAR FUNCTIONAL METHOD 
3.1 Introduction 

This chapter develops the main result of the dissertation. The functional method 
is a level of abstraction which permits complicated analyses of data with minimal 
effort. To make the presentation more straightforward, assume the data perfectly 
match the model 

/(*) = afce6fct> 

k - 1 

with the samples taken at integer times. This allows focussing on the problem as 
an interpolation problem even though more than 2 p data values may be used to 
determine 2 p unknown parameters. A later section presents methods which permit 
relaxing the specification so that the overdetermined inconsistent and the exactly 
determined interpolation cases can be solved. 

3.2 Theory 

Let $ be the set of all functions ip(t) : R — ► <C\0, of the form 

m = e w , 

which have the property that ip(t + 1) = 1). Let S be the space of all finite 

linear combinations of elements of <h. i.e. 

p 

mes => /(i) = £a*(l) (3.1) 

k = 1 

where pG Z is finite, G C, and 4>k{t ) E <F. 

Lemma 6. // Vi € S of the form in (3.1), and IF is a linear functional on S, then 

p 

Ui ) — ^ ak<f>k{i) [^4>k\ ( 0 ) 
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or, in matrix form with [P(f>k] (0) = \fF<j>k\ (^) U=o = A * 


F(y{i)) = [0i(O 02(0 ••• 0p(O] 

Proof. First note that 

Ui+ 1 = ai0l(O0l(l) + a 202(O02(l) + • • • + a p0p(O0p(l)> 

and because 0*(1) is constant, 

P{Vi) — a l [^01 ] (0 + a 2 02] (0 + 1- a p 0p] (0 



implies 

^■(?/i + i) = ai [JF0O (z)0i(l) + a 2 [J~ 02 ] (002(1) H (- a p [Pdp] (O0p(1)- 


So that if 


^(v<) = [ Wl] (0 [^02] (0 



then 

^(jfi+i) = [ [^0i] (0 [^02] (0 • • • [^0p] (0 



02(1) 
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Now, induct from i = 0 



so that for any i, from the above, 



Lemma 7. The matrix 

01 (0) 02(0) ... 0p(O) 

01 (1) 02(1) ••• 0p(l) 

01 (p 1) 02 (P 1) ••• 0 P (P-!) 
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is a Vandermonde matrix and can be represented as 


1 1 ... 1 

Z\ Z2 * • * Zp 




or 

1 1 ... 1 

e bl e h ... e 6p 

gbi(p-l) e b 2 (p-l) _ e bp(p-l ) 

for some Zk E C and bk E C. Then , provided z k 7 ^ z ji when k 7 ^ j, the matrix is 
nonsingular. 


Lemma 8. If b\,b 2 , ■ ■ ■ ,b v are distinct, all ai,a 2 , . . . ,a p are nonzero, and if 


Vk = ^2 a j 

j = 1 


e bjk 


for k — 1,2, ... , then the matrix 


H 0 = 


Vo 

2/i ■ • 

2/p — 1 

2/i 

2/2 •• 


2/p — 1 

2/p •• 

• V2p-2 


is nonsingular. 


Proof. Let 


Vk - 


Vk 

2/fc+l 


1^2/fc+p-l 
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Then, the companion matrix G which relates y k = Gy k _ l , is nonsingular. Note that 


y k = Gy k _ x 


implies 


Vk = G k y 0 . 


Now, if H 0 is singular, one of the columns can be written as a combination of the 
remaining columns. Assume that 

p- 1 

Vk= a i y r 

j=0,j^k 

Then 

p-i 

G k y 0 = E a ’ G ’y° 

j=0JjLk 

which, with a k = —1, gives 

v - 1 

0 = Y^ a jG j - 

3 = 0 

This last result is a polynomial in G of order less than p which is equal to zero. 
According to the Cayley-Hamilton theorem, zero is one of the roots of the character- 
istic polynomial for G. But since the eigenvalues of G are {e b \e b2 , . . . , e bp }, none of 
its eigenvalues are zero. This is a contradiction. Therefore, H 0 is nonsingular. ■ 


Theorem 9. Given a set of data y 0 , y u . . . , y n , ■ ■ • . generated at integer times by a sum 
of p distinct elements in S , then the matrix 



Vo 

2/i • • 

• 2/p-i 

2/i 

2/2 -• 

2/p 

2/p— l 

2/p • • 

• 2/2 P -i 
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is nonsingular, and the matrix G, given by 


Hvo) 

Hvi) ■■ 

F{y v - 1 ) 

?(Vi) 

Hvi) •• 

• Hy r ) 

0 

Hv r ) ■■ 

■ Hvip-i) 


has the the eigenvalues {^(0i(O)), JF(0 2 (O)), . . . , ^^(O))}. 


Proof. Consider the data vector 


y 0 = 


2/o 

2/i 


2/p— 1 


From the definition of the function which generated the data y 0 , there is a vector 

a i 

ao 

a = , 

Op 

so that 


01 (0) 02(0) 

0l(l) 02(1) 

J 


ax 

a 2 


2/o 

2/i 

_0lG?-l) 02 (,P 1) •• 

1 

t-H 

1 


dp 


_2/p— 1_ 


V'a = y 0 . (3.4) 


Then since V is nonsingular, 


a = V' l y 0 . 
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Also, from Lemma 6, if 


Huo) 
Hu i) 

Hup- 0 


then 


V 


HM<>)) 


Hh( o)) 


^ P (o)) 


o = Vo 


or 


V' 


Ai 


A? 


A r 


a = Vo- 


Where, as before: X k = ^(^(0)). Then substituting in the “Prony Step,” 

Ai 


V 


A2 


V l Vo = i/0- 


Now, designate 



2/t 


HVx) 

Vx = 

2/i+i 

, and = 

F{Ui+ 1 ) 


2/i+p-l 


^"( 2 /t+p-i) 


(3.5) 


(3-6) 


(3.7) 
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Claim: (3.7) holds for all i. 

Let 

*i(l) 

= 02 ( 1 ) 

M i) 


and 



Then, note 

<px{i) 02(0 ••• 0p(O 

<pi(i + 1) 02^ + 1) 0p(® d” 1) 

0i(i+p— 1) 02(^+P — l) 0p(*+p— 1) 


so that, 


VT*a = 


and 

VT*Ao = y t . 


Then perform the “Prony Step” again: 

VrAr~ l V~ l y t = y t - 


But since P and A are diagonal, 

TAf- 1 = ATT' 1 = A, 


= VTA 
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so 


VAV l y, = y x 


Let 

G = VAV' -1 . 


Then, 


Gy { = y l 


for each i € {0, 1 P - U Now, G takes each column of 

yo Vi ■ ■ ■ Vp-i 

yi V2 Vp 
Ho = . . . 

y P - 1 y P ■■■ V2 P -i_ 


to the corresponding column of 


Hvo) 

Hvi) 

. Hv P - 1 ) 

Hvi) 

Hvi) ■■ 

. Hyp) 

F(yp- 0 

Hv r ) " 

. F{y2 P -i) 


so 


GH 0 = H V 

Since G ~ A, the eigenvalues of G are the values 


{^(<*r(0)), ^(^(O)) ^(c>p(0))>- 
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Corollary 10. Given a set of data y Q , y x , . . . , y n , ■ ■ ■ , generated at integer times by a 
sum of p distinct elements in S, the matrix 

Vo 2/i • • • Up - 1 
2/1 2/2 ••• 2 /p 

2/p— i Up • • * 2/2p-i 

and the matrix 

r(yo) F(yi) ■ ■ ■ r(y p -i) 

Hvi) Hvi) Hvp) 

Hi — . . 

F{y P - 1 ) r(y p ) ... X{y 2 P -i) 

The generalized eigenvalue problem 

(Ho A - Hi)x = 0 

/ias the eigenvalues {iF((f>i(0)), lf(<j> 2 ( 0 )), ■ ■ • , H((t> p (0))}. 
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CHAPTER 4 
IMPLEMENTATION 


The previous chapter developed a functional theory which may be useful in find- 
ing the coefficients of the model equation. This chapter expands on this theme by 
developing some methods for using the functional approach to determine the coeffi- 
cients bk of the model equation. The first section describes the three basic categories 
of problem: the exact problem in which data are assumed to be perfect, the interpo- 
lation problem in which p data values and p functional values are used to determine 
an interpolating solution, and the overdetermined case in which there is an overabun- 
dance of data and functional values. The second section describes various types of 
functionals which have been identified to date. The third section connects some of 
the described functionals to methods which have appear elsewhere. 

4.1 Eigenvalues 

In this section, three different data situations are considered and approaches to 
using functionals to solve them are presented. The first situation is the theoretical 
one. In that situation, the data are assumed to be generated from the model equation 
and the data values are assumed perfect. Thus, it is called the “Perfect Data” case. 
Even in the case of perfect data, though, the numerical methods used to determine 
the coefficients, bk, may introduce anomalies, so two different methods of using the 
data to determine the coefficients are given. 

The second situation is the case in which there is “just enough” data to solve the 
problem. The developments in this case follow the method presented in Martin et al. 
[65]. The difficulty presented in this case is that, for some functionals, the matrices 
Hq and H\ are not full. Some means is needed to fill in the missing information 
to provide a complete solution to the problem. The results of this section do not 
have theoretical backing, however, the numerical results presented in Martin’s article 
indicate promise. 
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The overdetermined case is presented third. Standard linear least squares ideas are 
used to develop methods for coping with excess data. The theoretical foundations for 
using linear least squares for this type of approach are shaky however, and numerical 
results presented in a later section urge caution. 

4.1.1 Perfect Data 

Let x be an eigenvector of G J . Then, consider the equation 

H 0 G J = Hi, 

where H 0 and H x are two Hankel matrices related by a linear functional 
the right to obtain 

H 0 G t x = Hix 
H 0 Xx = H x x 

or 

(H 0 X - Hi) x = 0. ( 4 - 2 ) 

Thus, there are two methods of determining the eigenvalues of the matrix G. 
Method 1. Solve 

GH 0 = Hi 


. Apply x on 


(4.1) 


for G and determine the eigenvalues. 


Method 2. Solve the generalized eigenvalue equation 


{H 0 A - Hi)x = 0. 
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4.1.2 Interpolation 

Consider the situation in which there are p data values, y k , and p functional values, 
F(y k ) = y k . Then, to form the matrices and p - 1 unknown data values and 
p - 1 unknown functional values must be determined. The desired matrices are then 


y o 

2/i •• 

Vp-2 

2/p— i 

y i 

2/2 • • 

■ y v ~ i 

* 

y 2 

2/3 

* 

* 

2/p — 2 

2/p— i • • 

* 

* 

2/p — i 

* 

* 

* ., 

y o 

2/i • • 

• y P -2 

2/p-i 

2/i 

h ■ ■ 

. . 2/p-i 

* 

h 

2/3 • 

. . * 

* 

Vv- 2 

2/p— i • 

* 

* 

2/p— i 

* 

* 

* 


where the * represent unknown values. Then the solutions, A, to the eigenproblem 

(. H 0 \ - Hi)x = 0, 

are the values *(*(0)), “ usua1 ' The difficulty UeS in determining Unkn0KD 
parts of the matrices. 

It may be possible to create an iterative method which will converge to the correct 
matrices H u i — 0,1. Suppose for a moment that the values y p ,y P +u ■ ■ ■ , Vip-2 ancl 
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assumed to be zero. Then let 


2/pt Vp+l, • • • > 2/2p-2 ^e 



y o 

2/i • • 

• y P -2 

Vp- 1 


y i 

2/2 •• 

• y P -i 

0 

o 

II 

y 2 

ya •• 

. 0 

0 


Vp- 2 

2/p-i • 

0 

0 


Vp— i 

0 . 

. 0 

0 


y o 

2/1 • 

■ • Vp- 2 

2/p-i 


V i 

2/2 • 

■ ■ y P -i 

0 

Jjf> = 

2/2 

ya • 

.. 0 

0 


Vp- 2 

2/p-i 

.. 0 

0 


2/p-i 

0 . 

.. 0 

0 


and solve the eigenproblem 

{H 0 X - Hi)x = 0, 

for the values A^. 0) . Solve the functional equation 

HMO)) = 4 0) 


for 4> k (l) = z k and form the polynomial 

(<°H*) = no - 

Jfc=l j= o 


and form the companion matrix 


4(0) = 


0 

0 


1 0 
0 1 


— G () — &1 


0 

0 


dp— i 
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Then, since both y k and y k are sequences which obey the same difference equation 1 , 
use the matrix A to fill in a first approximation to the missing values of H 0 and 
Hl . Call the new matrices H and h[ 1] . Repeat the process until convergence is 

achieved. 

This method has been presented without proof. However, in experiments con- 
ducted by Martin et al. [65], the method converged in only 2 or 3 iterations for the 
oscillatory (derivative) functional. An example of this method in use is shown m 

section 5.3. 

4.1.3 Approximation 

The formulation of the functional approach assumed that there was exactly enough 
data to create square matrices which could be inverted. As it happens, if there is 
an excess of data, and the matrix G represents a valid relationship between the data 
and the transformed data, then its eigenvalues still represent the functional values. 
Then the only concern is in determining the eigenvalues of G in a mathematically 

valid manner. 

However, in the case where the system 

GHl = Hi 

is overdetermined and inconsistent, some means must be found which offers a reason- 
able approximation of the eigenvalues of the problem. The most obvious method of 
determining the eigenvalues is to determine the G which minimizes the norm of 

GHl - Hj , 

and then to determine the eigenvalues of G. One way to do this is to use linear least 
squares to determine the rows of G. The use of the normal equations gives 

Hj H 0 G J = Hq Hi. 

x The linear functionals pass across the difference equations. 
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The solution, G, is unique if the matrix Hj H 0 is nonsingular. 

The numerical results reported in a later section, page 37, suggest that this for- 
mulation of the problem, for any method of solution, is inconsistent. A much better 
formulation of the problem is to use the generalized eigenvalue problem. 

Method 3. Let 


y o 

y i 

to 

r- 

i-H 

1 

Vi 

y 2 

2/3 

2/p 

Z/n— p— r+ 1 

2/n— p— r+2 

2/n-p-r+3 

3 

1 

1 

o 

Vi 

2/2 

V P - 1 

y i 

Vi 

2/3 

yp 

jjn-p-r+l 

y-n—p—r+2 

2/n— p— r+3 ■ ■ 

Vn—p—r 


Then solve 

(HjH 0 \-HjH l )x = Q 

as a general eigenvalue problem. 

4.2 Functionals 

In this section, three different types of functionals are presented. The shift func- 
tionals are seen to be variations on the shift operator. The discussion shows that 
there are many different ways to use shifts to obtain solvable linear systems. The 
derivative operator is also a linear functional and will be seen to offer a simple way 
of determining the coefficients. Finally, the anti-derivative or integral operator is 
another linear functional which offers some simple solution methods. 
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4.2.1 Shift Functionals 

Consider the shift operator qy(t ) = y(t + 1). Then, it is a small step to create a 
linear functional from the shift functional by creating a polynomial of shifts, 

T{y) = (c 0 + Ciq + c 2 q 2 + • • • + c T q T )y. 


Let 

y n = T{y n ) — (c 0 + c x q + c 2 q 2 -I 1- c T q r )y n - 

Form H 0 from the values y n , and Hi from the values y n . Then the solutions, 
Ai, A 2 , . • • , A p , to 

{Ho A - Hi)x = 0 ( 4 - 3 ) 


are such that 

A fc = c 0 + cie 6fe + c 2 e bk2 + • • • + c r e hkT 
Afc = Cq + C\Zk + C 2 ^ + • • • + C r z k 
0 = Co - Afc + c x z k + c 2 z\ + • - • + CrZ T k , (4.4) 

where b k is one of the true coefficients. The polynomial (4.4) has r roots, of which 
there is only a guarantee that at least one is the desired root z k . Thus there is some 
ambiguity in choosing the coefficient b k based on the polynomial (4.4). However, there 
is some hope. Under the assumption that the coefficients b k are negative real, then 
the desired root, z k , must be in the interval (0, 1). Moreover, the following lemma 
and corollary show that if A obeys certain easily tested criteria, and the polynomial 
is of a certain form, then it is easy to see that there is exactly one root of (4.4) in 

( 0 , 1 ). 

Lemma 11. If0<a<r,the polynomial equation 

a — z + z 2 4 1 - z r 

has exactly one solution in (0,1). 
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Proof. Let 


p(z) = Z + Z 2 + • • • + Z T 


then. p( 0) = 0, p(l) = r and z > 0 =► p'{z) > 0. 
Therefore, p(z) : (0, 1) - (0, r) is one-to-one and onto. 

Corollary 12. // 1 < a < r + 1, 


OC = l + Z + Z 2 -\ 1 - Z T 


has exactly one solution in (0, 1). 

Lemma 11 can be extended to polynomials with all positive coefficients if desired. 
This alleviates some of the difficulties associated with the use of shift polynomials for 
determining coefficients, but requires certain a priori assumptions to be made about 

the nature of the coefficients. 

There are certain major advantages of using the shift functionals. Consider a 
system with coefficients 6, = -0.0202 and b 2 = -0.0101 which correspond to 0.98 = 
g— 0.0202 0.99 = e -0 0101 respectively. Then, if the coefficients are to be determined 

from data which is generated by the function 

/(t)=c 1 e-° 0202t + c 2 e-°- olou , 

the shift functional T = 1 + 9 + T + • • • + provides a characteristic equation with 
roots 6.5937 and 6.7935. The solution of this equation is more numerically stable than 
the equation with roots 0.98 and 0.99. Since the objective is to solve the generalized 

eigenvalue problem 

{Ho A - Hi)x = 0 

for the roots, then the use of the functional to add separation to the roots is beneficial. 
The tradeoff is that the polynomial 

0 = c 0 - A fc + ciz fc + c 2 z 2 k + h c r z T k 
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is hard to solve. If this is a large consideration and enough data is available, then 
instead of polynomials of shift, the user can just use a shift of high degree. To see an 
example of one of these shift functionals in use see page 26. The functional approach 
allows for easy justification of the resulting method and easy calculation of the inverse 
function. For comparison with the above, the functional transformation T(y) — <J y 
yields corresponding eigenvalues of .784/ and .8864. 


4.2.2 Differentiation Functionals 


Let 


™ = Jt v ■ 


Then, the eigenvalues of the matrix G are 


HMO)) = !«*“ 


= b k . 

t = o 


Then, provided the derivatives of the measured data are available, this formulation 
presents a very simple method of determining the decay factors. The theoretical 
groundwork laid in developing the functional technique, significantly reduces the al- 
gebraic overhead in describing the use of derivatives. Compare this to the development 
in [65] of osculatory interpolation. 

Often the determination of derivatives from data is difficult. Noise in the data 
can significantly affect the computed derivatives, it is easier to compute the integral 
of the data. The functional technique can be used to justify this approach also. 


4.2.3 Integration Functionals 



The second integral is in place to allow the integrated process to be shifted to satisfy 
the stationary assumption of the AR process. Assume that Re b k < 0. Then the 
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eigenvalues of G are 


/ oo 


e bkT dr 


t=0 



£ 

S' 


One possible advantage to using the integration functional is that integration is 
an inherently smoothing operation and there may be some beneficial effect on the 
noise content of the data by using numerical integration. 

4.3 Comparison with Methods from Literature 
The functional technique can also be viewed as a unifying theory for various 
methods that exist in the literature. Prony’s method and the matrix pencil are 
the most obvious applications of the functional method. The osculatory method 
developed by Martin et al. is another method. 

Let q be the standard shift of degree one. Then put 

?{v) = qy- 

Let Ho and H x be as usual. Solve the matrix equation 

GHo — H x . ( 4 - 5 ) 

Then, the matrix, G, is a companion matrix whose last row, [-c 0 , -c u . • • , -c p ], is 
precisely the vector which results from Prony’s method. The eigenvalues of G are the 
eigenvalues resulting from the matrix pencil and are the roots of the characteristic 
equation identified by Prony’s method. Now, consider the problem as an eigenprob- 
lem. The eigenvalues of G are theoretically the same as those from 

(HqX — H x )x = 0, x ^ 0- 

The only time that a difference might occur in the solutions to the two problems 
is when the conditioning of the matrix H 0 is such that (4.5) is unstable [30, 106], 
Therefore, the simple first order shift, F{y) = W. generates Prony’s Method and the 
matrix pencil. Examples of these methods can be found in [2, 35, 63, 79]. 
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Let T(y) — j t V- Designate the derivative values by y n . Then let H 0 be given 
as usual for data values y 0 ,y u . . • , y 2p - 2 and let H x be generated from the derivative 
values jfe, Vu • • • , 2/2p-i- Then the solutions, A, to the eigenproblem 

(Hq\ — Hi)x = 0, x ^ 0, 

are precisely the coeffcients, b k . This result is the same as that achieved by Martin 
et al. [65]. 

Thus, the standard AR formulation of the problem is the special case of the func- 
tional method resulting from the shift functional of order one. The new development 
of Martin et al. is the special case of the functional method resulting form the deriva- 
tive functional. 


4.4 Conclusion 

The implementation of the functional method is simple. It only requires that a 
linear functional be chosen for which the inverse problem can be solved. Then, the 
data is transformed using the functional and the eigenvalue problem solved. Then 
the user solves the inverse linear functional problem for the coefficients of the system. 
The use of the derivatives in one case showed that if there is data available which is 
related to the measurements by some linear functional, that additional data can be 

used to aid the process of determining the coefficients. 

This chapter showed some simple functionals which are easy to implement with 
measured data. In the case of the shift-polynomials, the inverse problem is not neces- 
sarily easy to solve, however, the division of difficulty between the eigenvalue problem 
and inverse functional problem may prove to be beneficial in certain instances. 

The integration technique has promise in the area of reducing the sensitivity of 
the method to noise. The comparisons in the next chapter between methods based 
on different functionals shows this promise clearly. 
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CHAPTER 5 
NUMERICAL RESULTS 


The last chapter demonstrated some simple linear functionals which can be used 
as bases for solution methods. This chapter compares some of those functional based 
extensions of Prony’s method to solution methods which are reported in literature. 
The results reported here show that the functional methods are not uniformly better 
than methods which already exist. However, some of the methods presented here 
show promise in certain situations. The lesson to be learned from these examples is 
that the functional approach offers the user the ability to try different methods until 
he is satisfied with the results. The addition of many new methods can only result in 
better analysis capability. 


5.1 A Simulation 

This section presents the results a simulation comparison of several AR meth- 
ods versus several different methods based on linear functionals. First, consider the 
problem of identifying the nonlinear coefficients of the function 

f(t) = e~ l + e~ 2t + e~ 3t + e~ 4t + e~ 5t 

from data obtained by simulating the function at the evenly spaced times, U — ih for 
i 6 {0, 1, ... , 100} and h = 0.01. 

Several different methods are used to try to solve this problem. Pisarenko Har- 
monic Decomposition (PHD), standard least squares, and Osborne’s Objective Re- 
weighting Algorithm (ORA) are representative algorithms from the AR class of tech- 
niques. Several different methods based on the functional approach are used for the 
remainder of the comparison. For the Derivative functional, the actual function was 
used to generate derivative values with no added noise. For the integration functional 
used in the two examples “Integrate-1” and “Integrate-2” , the eigenvalue problems 
were HjH 0 A - H^H U and HjH x A - Hj H 0 respectively. The “Shift-n” functionals 
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are the simple shift of order n functionals. These provide easy inverse solution with 
no am biguity in the solution. All of the functional methods were solved as eigenvalue 
problems with the QZ algorithm. 

Table 5.1 shows the results of the simulation when the level of noise in the data 
is strictly a function of the numerical noise in the computer representation. For 
Table 5.2, measurement noise was simulated in the data at the level of N(0, .00005 2 ). 
Certain data values represent a negative root to the AR equation. Some of these values 
are easy to spot and are designated with *. Note that for the Shift-18 functional, 
17.5 x 18 = 315, which is 7 r/h within rounding. The same is true for the Shift-21 
functional, where 15 x 21 = 315. 

5.2 Lanczos’ Problem 

Lanczos presented an example of using Prony’s method [56]. In his example, he 
summed the data to create a smaller data set which could be solved directly using 
Prony’s method. The data was generated by the function 

f(t) = 0.0951e -t + 0.8607e _3t + 1.5576e“ 5 * 

with equal length time intervals h = 0.05. Each data value was rounded to 2 decimal 
places to simulate measurement capabilities. The data are reproduced in Table 5.3 
for reference. 

Both Lanczos and Ruhe [89] report limited success determining the exponential 
coefficients. Table 5.4 reports the results obtained by Lanczos, Ruhe, and current ex- 
periments using Osborne’s Objective- function Reweighting Algorithm(ORA) method, 
Linear least Squares Prony method, the Matrix Pencil, and several different functional 
methods. The factors, k, represent the model size used in the solution example. Note 
that in that table, several of the examples use oversize models. If the user is interested 
in only real decay factors, the use of an oversize model can increase accuracy in the 
real factors. This is most apparent in the results reported for k = 11. Outside of 
the ORA method of Osborne, this functional returned the best results for this data 
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set. In the case of the linear least squares and matrix pencil methods, some of the 
returned coefficients were complex. These discarded answers are represented in the 
table by a dash. 

The method developed by Osborne and reported in [42] as the ORA method is a 
simplification of the Gradient-function Reweighting Algorithm (GRA) which was not 
used for these comparisons. The ORA method appears to provide answers which 
are significantly better than those available from any other method. These results 
can be misleading however. ORA is a nonlinear method which requires a starting 
value for the autoregressive parameters. Global convergence is not guaranteed, so 
it is customary to start the algorithm with the true parameters. Alas, if the true 
parameters are not known, the algorithm may not converge to the correct values. 

The result from the functional method T\ with model size k=ll is also good. For 
reference, the functionals used in Table 5.4 are listed here. 

T\ = 1 + 2 . 

JF 2 = 1 + z + z 2 + 2 3 + z 4 . 


The functional T\ is invertible, but is not. For any given solution, A, of the 
eigenproblem, the polynomial equation 

A = 1 + z + z“ + z 3 + z 4 

has 4 solutions. The difficulty of the method is that it must distinguish among these 
solutions and report the one which is most likely. The heuristic used in this program 
is that the correct solution is the one which is closest to one on the complex plane. 
Imaginary parts are then discarded and duplicates merged. The zeros reported in the 
table represent failures to find any solution for at least one A which was within 1 unit 
distance from 1 on the complex plane. 

At then end of the table, results are reported for several different integration 
methods. All of the integration methods use the same type of integration. That is, 
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they add the values in sequence, generating a new data stream of cumulative sums. 
The entire sequence then has the last data value subtracted from it to simulate an 
integration from x to infinity (see page 33). This works well, as can be seen from 
the results. The different results stem from the method used to solve the eigenvalue 
problem. The problem is overdetermined so the projection method shown in section 
4.1.3 is used. In the first two examples, Integrate-1 and Integrate-2, the eigenvalue 
problem is 

Hj HqG = Hj Hi Integrate-1 


or 


{HjH 0 X - HjHi)x = 0 Integrate-2. 

Notice that this differs from the usual projection in that the matrix from the modified 
data is used to define the projection space. The results from both solution methods 
are shown to be the same to within reported significance. The next two examples are 
made from the usual projection. 


Hi H 0 G = Hi H x Integrate-3 

and 

(Hi H 0 A - Hi H x )x = 0 x^O Integrate-4. 

In the case reported as Integrate-3, the equation is solved for G and the eigenvalues 
are found from G. In Integrate-4, the generalized eigenvalue problem is solved directly 
using the QZ algorithm described in [30]. Clearly, using H\ to define the projection 
space greatly improves the conditioning of the problem. One possible explanation 
for this is that the process of summing the data reduces the influence of rounding 
errors. The matrix HjHi is thus better conditioned than Hj H 0 . However, the 
results for Integrate-3 and Integrate-4 clearly show that solving the eigenproblem 
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directly is much more stable than solving the matrix equation first. In the event that 
projection onto the space defined by Hi is not justified (all the shift functionals), the 
eigenproblem should be solved directly to avoid numerical difficulties. The user may 
even want to experiment with the more expensive VZ algorithm reported in [106]. 


5.3 An Interpolation Example 

Suppose data were collected from an experiment and due to some difficulty, the 
middle part of the data were deemed inadequate. For the sake of definiteness, suppose 
the data values y 3 , 7/4 of a total of 8 data values were missing. Then, the researcher 
has only t/ 0 > Vi, V 2 , 2 / 5 , 2/6, 2/7 from which to make a determination of the decay rate of 
the experiment. Since the data can be taken in pairs, 2/fc,2/Jfc+5, the researcher chooses 
to use a shift functional of order 5 to determine the coefficients, in this case, the 
starting matrices look like 

Vo 2/i 2/2 

Ho = 2/1 V2 0 

V2 0 0 _ 

and 

2/5 2/6 2/7 

Hi = y 6 y 7 0 . 

.2/7 0 0 _ 

To see how this example actually performs, test data were generated with 

f{t ) = e -t + e~ 3t 4 - e~ 5t , 


for the interval [ 0 , 1 ] at equal length time intervals h — 0.05. No noise was added. 
Table 5.5 contains the first ten data values to 4 decimal places. It also has the results 
of the interpolation after 6 iterations of the scheme. The data values used in the 
experiment were maintained in their full precision in Matlab on a PC. After three 
iterations the matrices are 

’3.0000 2.5907 2.2522 
H 0 = 2.5907 2.2522 1.9717 
2.2522 1.9709 1.7381 
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and 


H 1 = 


1.5377 

1.3705 

1.2284 


1.3705 

1.2284 

1.1068 


1.2284 

1.1070 

1.0026 


and the coeficients, b k , have been estimated as {-1.2609,-4.5082,-6.4569}. Note 
that after three iterations, the matrices are not true Hankel matrices and the proper 
values for the interpolated data cannot be determined. After three more iterations, 
the matrices are truly Hankel matrices to 4 decimal places and the estimated data 
stream is reproduced in the second column of Table 5.5. The values with * be- 
side them represent the interpolated values. The estimated coefficients are now 
{—1.2638, —3.1055, —4.5049}. Thus the method appears to converge nicely and in- 
terpolates the missing data values to within 0.0003. Further iterations result in little 
to no improvement. The Matlab program which implements the iteration is listed in 
Figure 5.1. 

This example shows that the functional method can be useful in interpolating 
data. Caution must be used; however, the same procedure was tried with the first 
ten values from Table 5.3. The iteration scheme never converged to Hankel matrices 
with 4 decimal places of accuracy (to be expected from data with only 2 decimal places 
of accuracy), although the returned coefficients were closer to the correct values than 
many of the examples in Table 5.4. 

A note about implementation is in order. The algorithm in Figure 5.1 specifies 
that the columns in the matrices Ho and H\ be updated from the preceding columns 
in the previous iteration. In other words, if r is the k - 1st column in Hq\ then the 
fcth column in Hq + ^ should be 


G (l) r 

instead of 


(5.1) 


(G^-'y 


(5.2) 
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where y is the fully known first column of Ho- The update method (5.2) results in 
all of the coefficients converging to the same number rather than (in this case) three 
distinct ones. To gain some insight into why this is so, consider a step of the iteration 
using (5.2): 




V (G<‘>) 2 y 


(5.3) 


and a step of the iteration using (5.1): 

= [y G {i) y G^G^y] (5.4) 


It appears that slowing down the convergence rate by using a previous value for 
G keeps the algorithm from overconverging. It should also be noted that with the 
matrix G in companion form, the known parts of H 0 and Hi never change. 
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Table 5.1: Results with no added noise 


Method 

Results 

PHD 

56 + 230 i 

56 - 2302 

-1.09 

-2.77 

-4.76 

ORA 





-4.73 

Derivative 

-5.0000 

-4.0003 

-3.0005 

-2.0001 

- 1.0000 

Shift- 1 

-5.0016 

-4.0161 

-3.0223 

-2.0060 

-1.0002 

Stand. LLS 

-5.0016 

-4.0163 

-3.0226 

-2.0060 

-1.0002 

Shift-3 

-4.9995 

-3.9952 

-2.9929 

-1.9981 

-.9999 


Table 5.2: Results with noise level e ~ iV(0, 0.00005 2 ). 


Method 

Results 

PHD 

-1.56 

-5.31 

-6.50 

-4.5 + 1672 

-4.5 - 1672 

ORA 

-1.08 

-4.71 

-2.71 

-5.98 - 3082 

111 - 102 

Derivative 

-1.43 

-4.21 

-0.28 

0.0000 

0.0000 

Integration-1 

-1.09 

-5.76 

-2.25 

-173 - 682 

-173 + 682 

Integration-2 

-1.59 

-4.26 

0.23 

-142 + 522 

-142 - 52i 

Shift- 1 

-1.37 

-4.18 

-60 + 314i* 

-47 + 1842 

-47 - 184i 

Stand. LLS 

-1.37 

-4.18 

-60 + 314i* 

-47 - 1842 

-47 + 1842 

Shift-3 

-1.37 

-4.17 

-102 

-36 - 562 

—36 + 56i 

shift-6 

-1.36 

-4.15 

-40.8 

-23.0 + 26i 

-23.0 - 26z 

shift-9 

-1.36 

-4.14 

-19.9 

-19.0 - 34.92 

-44.1 - 34.9i 

shift- 12 

-1.35 

-4.13 

-32.6 

— 22.0 ~b 9.0/i 

-22.0 - 9.072 

shift-15 

-1.35 

-4.11 

-24.6 

-14.7- 3.17i 

-14.7 + 3.172 

shift- 18 

-1.35 

-4.11 

-15.5 + 17.52* 

-8.89 + 16.4i 

-8.89 - 16.4i 

shift-21 

-1.34 

-4.09 

— ( .58 + 15.02* 

-6.81 + 6.032 

-6.81 - 6.032 

shift-24 

-1.34 

-4.06 

-5.44 

-3.66 + 2.942 

-3.66 - 2.94i 
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Table 5.3: Lanczos’ data. 


Data Values 

2.51 

1.12 

0.53 

0.27 

0.15 

0.09 

2.04 

0.93 

0.45 

0.23 

0.13 

0.08 

1.67 

0.77 

0.38 

0.20 

0.11 

0.07 

1.37 

0.64 

0.32 

0.17 

0.10 

0.06 


Table 5.4: Results from various methods. 


Method 

Result 

Ruhe (k=2) 

-1.75 

-4.55 


Lanczos (k=2) 

-1.58 

-4.45 


mmmm 

-3.5 

-8.5 


LLS (k=3) 

-2.71 

-5.1 

— 

LLS (k=4) 

-1.97 

-4.6 

— 

Matrix Pencil (k=2) 

-3.53 

-8.5 


Matrix Pencil (k=3) 

-2.71 

Oi 

o 

OO 

— 

ORA (k=2) 

-1.81 



ORA (k=3) 

-0.97 


-5.1 

PHD (k=3) 

-0.26 

-3.92 

-16.9 

(k=3) 

-2.71 

-5.08 

0 

^ (k=ll) 

-0.89 

-3.35 

-5.29 

(k=2) 

-2.8 

-5.3 


(k=3) 

-2.0 

-4.7 

0 

Integrate- 1 (k=3) 

-.72 

-3.2 

-6.5 

Integrate-2 (k=3) 

-.72 

-3.2 

-6.5 

1 Integrate-3 (k=3) 

-.3.3 

-16.2 

-65.5 


-.3 

-1.2 

-6.0 
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Table 5.5: Comparison of true data with interpolated data 


Actual 

Interpolated 

3.0000 

3.0000 

2.5907 

2.5907 

2.2522 

2.2522 

1.9707 

1.9708* 

1.7354 

1.7357* 

1.5377 

1.5377 

1.3705 

1.3705 

1.2284 

1.2284 

1.1068 

1.1068* 

1.0023 

1.0020* 


7, assumes the polynomial p and the matrices HO, HI already exist 
G=[0 1 0 
0 0 1 

-p(4 : -1 : 2)] ; 

7, it is very important to update HO, HI from HO, HI 
7. rather than from y, yl and G 
H0= [y G*H0( : , 1) G*H0(: ,2)] 

Hl= [yl G*H1 ( : , 1) G*H1(: ,2)] 
e=real(eig(Hl ,H0) ) ; 

el=log(e)/5; 7. assumes yl is shifted 5 from y 
p=poly(exp(el)) ; 

e2=el/.05 7. only necessary when reporting coefficients 

Figure 5.1: Matlab Program for interpolation 
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CHAPTER 6 
CONCLUSIONS 


The functional technique described in this dissertation represents a great leap 
forward for the task of solving Prony’s Problem. Before this technique was developed, 
the available methods depended on either nonlinear techniques or the autoregressive 
formulation. The AR formulation has been the method of choice because the problem 
is linear and responds to many different techniques. With the introduction of linear 
functionals, the arena has become much more vast. 

This dissertation showed the theoretical basis for the functional techniques. It 
also showed how simple many problems become when the functional techniques are 
used. The justification for many different approaches to the problem becomes easier 
and the different transformations have been unified with one theory. 

However, this dissertation could only scratch the surface of the new functional 
technique. There are still many questions to be answered. The iteration technique 
introduced in section 4.1.2 needs justification. The limitations on its convergence 
need to be identified. Since those convergence properties depend on the functional 
which is used, it is impossible to plumb those depths in this dissertation. Addition- 
ally, the properties of the overdetermined case need to be investigated. The method 
presented there is apparently irreconcilable with the covariance methods presented 
by Pisarenko, Osborne and others [83, 73], Hopefully further research in this area 
will result using the functionals to create an orthogonalization technique for these 
problems. The existence of such a technique would result in vastly improved compu- 
tational performance in these problems. 
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